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Abstract 

Systems of differential-algebraic equations (DAEs) represent a widespread formalism in 
the modeling of constrained mechanical systems and electrical networks. Due to the 
automatic, object-oriented generation of the equations of motion and the resulting re¬ 
dundancies in the descriptor variables, DAE systems often reach a very high order. This 
motivates the use of model order reduction (MOR) techniques that capture the relevant 
input-output dynamics in a reduced model of much smaller order, while satisfying the 
constraints and preserving fundamental properties. Due to their particular structure, 
new MOR techniques designed to work directly on the DAE are required that reduce the 
dynamical part while preserving the algebraic. In this contribution, we exploit the spe¬ 
cific structure of index-1 systems in semi-explicit form and present two different methods 
for stability-preserving MOR of DAEs. The hrst technique preserves strictly dissipativity 
of the underlying dynamics, the second takes advantage of 'H 2 -pseudo-optimal reduction 
and further allows for an adaptive selection of reduction parameters such as reduced 
order and Krylov shifts. 
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1. Introduction 

The accurate modeling of dynamical systems such as flexible mechanical structures, 
electrical circuits with large-scale integration and interconnected power networks often 
results in mathematical models of very high dimension, meaning that the number of state 
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variables and equations used to represent the dynamic behavior grows to the order of 
tens of thousands up to millions. 

This curse of dimensionality is even more accentuated if the formalism used to model 
the system yields a set of differential-algebraic equations (DAEs). This type of formu¬ 
lation arises e.g. if the system is modeled by equating the fundamental physics of each 
individual element, which is then linked to its neighbors through appropriate intercon¬ 
nections. This approach has the fundamental advantage of being quite generic, making it 
therefore suited for automated, object-oriented computerized modeling. DAEs represent 
systems of equations that are composed of a dynamical part, describing how the state 
can vary in time, and an algebraic part, describing static relationships and constraints 
amongst state variables. Subsequently, within this formalism the states used to model 
the system do not represent minimal coordinates needed to describe the dynamics but 
generally include a set of algebraic variables used for expressing a constraint manifold on 
which the state of the system has to evolve at all times. Due to this additional algebraic 
variables, the number of the states used to model the system becomes even larger. 

Especially in a large-scale setting, it is not computationally easy—if at all possible— 
to separate the dynamic states from the algebraic ones. This implies that the DAE has 
to be analyzed and treated directly, motivating the flourishing theory on DAEs that has 
been developed in the past decades [1-3]. 

In the following contribution, we shall focus on the simpler, yet common case of linear, 
constant coefficients DAEs given as a generalized state-space system of the form 

E X = Ax B u 

( 1 ) 

y = C X D u 


where E G is the singular descriptor matrix, A G is the system matrix and 

xGR^, uGR'”, yGR^ {p,m<^N) represent the state, input and output of the system 
respectively. Whenever the matrix E is regular, we will refer to (1) as a system of 
ordinary differential equations (ODE) in implicit form. 

If the model in (1) is used for simulation, optimization or control synthesis, a high 
order N might pose severe limitations in terms of computation time and even storage 
capabilities. For this reason, model order reduction (MOR) techniques are implemented 
to obtain a significantly smaller reduced order model (ROM) that captures the dominant 
input-output behavior of the full order model (FOM) while preserving characteristic 
properties. 

A common approach to the reduction of state-space models like in (1) is given by the 
Petrov-Galerkin projection =Ii = EV {W^ E V) ^ and yields 



A,. B,. 

W^AV Xr + W^B u 
CW Xr -k U 

Cr Dr 


( 2 ) 


where Xr € M” (n <C N) represents the reduced state vector. Accordingly, the aim of 
MOR techniques lies entirely in the appropriate choice and design of projection matrices 
V,W. 
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MOR theory for linear state-space systems satisfying det E^O has been extensively 
studied over the past decades and includes established techniques such as modal reduction^ 
balanced truncation and moment matehing [4]. The latter, also known as Krylov-sub space 
methods, implicit moment matching or interpolatory reduction techniques, stand out due 
to their generality and low computational cost, making them a predestined candidate for 
the reduction of truly large-scale systems. However, these reduced computations reveal 
only local information about the system and global properties of the model might be 
lost. For instance, one challenge arising when applying moment matching is that, given 
a stable FOM, stability of the ROM is not guaranteed per se. While there are some 
approaches to preserve stability in the case of ODEs, stability preservation in the more 
general case of DAEs has—to the authors’ knowledge—not been addressed so far. 

In this contribution this problem is considered. Two different techniques for stability 
preservation in interpolatory reduction are presented. Both techniques apply to index-1 
DAEs in semi-explicit form and do not require the explicit computation of the underlying 
ordinary differential equation. Further it will be shown in theory and through numerical 
examples how an orthogonal projection (V = W) of such DAEs can be applied correctly 
only if some conditions are met. Einally, we will explain how to adaptively choose reduced 
order and interpolation points of the Krylov subspaces for this class of systems. 

The sequel of this contribution is outlined as follows: section 2 revises the existing 
theory on model reduction, stability preservation and DAEs. In section 3 a first stability 
preserving reduction technique is presented. This requires the extension of the concept of 
strictly dissipative realization of a state-space system to the more general case of DAEs 
and shows how to exploit this property for stability preservation. The correct reduction 
of DAEs by orthogonal projection will also be discussed. Section 4 presents the second 
stability preserving technique based on 'H 2 -pseudo-optimal reduction. For this purposes, 
the Pseudo-Optimal Rational Krylov (PORK) algorithm will be adapted. Further, the 
adaptive choice of reduced order and interpolation points will also be generalized. Finally, 
the theoretical results are validated through numerical examples in section 5. 

2. Preliminaries and problem statement 

2.1. Model reduction of ODEs by Krylov-sub space methods 

As introduced in section 1, the computation of a reduced model by means of pro¬ 
jection boils down to finding suitable projection matrices V,W (cf. equation (2)). De¬ 
pending on which properties of the original model we would like to preserve, several 
methods of determining the projection subspace have been studied in the literature. For 
instance, modal reduction is aimed at preserving dominant eigenmodes of the system [5], 
while balanced truncation preserves dominant directions in the system that are equally 
controllable and observable [4]. Due to their generality and low computational cost, 
Krylov-sub space methods have been studied intensively over the past decades [6-8] and 
shall be introduced briefly in the following. 

If the matrices V,W in (2) span n*^-order tangential input and/or output Krylov 
subspaces 


Kn {{A - soE)-^E, (A - soE)-^B r) 

ICn {{A - soE)-^E^, (A - soE)-^C^ l) 
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(3a) 

(3b) 



respectively for some complex shifts Sq G C and tangential directions r G K™, I G Mp, 
then the reduced model matches some coefficients of the Taylor series expansion of the 
original transfer function around sq- In the case of systems with more than one input or 
output (m,p>l) the Taylor coefficients are not scalar and the matching applies only to 
the linear combination defined by the tangential directions r and 1. The computation of 
such matrices requires only the solution of (typically sparse) linear systems of equations, 
which makes it particularly attractive in a large-scale setting. In this framework, the task 
of computing good reduced models translates to hnding suitable shifts and tangential 
directions. 

For theoretical considerations it is worth noting that, at least in the ODE case, there 
is a duality between the Krylov subspaces in (3) and the sparse-dense Sylvester equations 

AV -E VSv -BR = 0 (4a) 

W^A-SwW^E-LC = 0 (4b) 

respectively, where the pairs (Sy^R) and {Sw,L) encode the interpolation data (shifts 
and tangential directions) [9, 10]. 

While reduction methods like modal techniques and balanced truncation intrinsically 
preserve stability, this is generally not the case for Krylov-subspace methods. This can 
be justified by the fact that while the former methods require costly computations that 
reveal the structure of the system (eigenvectors or Gramians respectively), the latter 
method limits the effort to the computation of only a few directions that guarantee 
matching of the transfer behavior at certain frequencies. It is therefore highly desirable 
to include some strategy to ensure stability of the ROM even for interpolatory model 
reduction nonetheless. 

2.2. Stability preservation in MOR of ODEs 

A few results are available for the stability preservation in interpolatory reduction of 
ODEs and shall be revised in the following. 

2.2.1. Strictly-dissipative realizations 

A state-space system as in (1) satisfying detAT^O is said to be in a strictly dissipative 
form if and only if E = E^ ;^0 and AG-A^ ^0. Since the trajectory of such a system can 
be shown to be strictly decreasing with respect to the elliptic 2-norm induced by E [11], 
i.e. 

Mt)\\lE ■= x^{t)Ex{t), 

systems in strictly dissipative form are asymptotically stable. In other terms, E = E^ >- 0 
and A-\-A^ -< 0 implies that all generalized eigenvalues of the pair (A, E) have strictly 
negative real part. This fact can be exploited in a projective MOR setting as in (2) by 
performing orthogonal projection with a full-column-ranked matrix V = W. This type 
of (Galerkin) projection preserves dehniteness of the original matrices and it is therefore 
easy to show that the reduced model in (2) will preserve strict dissipativity and hence 
be stable [12]. Please note that any asymptotically stable ODE can be transformed to 
a strictly dissipative realization and that systems in second order structure can often be 
directly transformed into an equivalent strictly dissipative first order realization [11]. 
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2.2.2. 'H 2 -pseudo-optimality and CUREd SPARK 

Recently, interpolatory techniques for ODEs that adaptively choose the Krylov shifts 
and guarantee stability preservation at convergence have been proposed [11, 13, 14]. In 
this contribution, the focus will lie on the Stability-Preserving, Adaptive Rational Krylov 
(SPARK) Algorithm introduced by Panzer et al. in [14] and further developed in [11] 
due to its suitability in the reduction of DAEs. The stability preservation in SPARK 
is guaranteed by construction, exploiting the concept of 'H 2 -pseudo-optimality [10, 15]. 
Given the desired interpolation data encoded in (S'y,i?) (cf. equation (4)), it is possible 
to directly construct a reduced model which satisfies the interpolatory conditions and 
shares the same spectrum as {—Sy). Based on the Sylvester equation (4), the Pseudo- 
Optimal Rational Krylov (PORK) algorithm yields the desired pseudo-optimal reduced 
model [15]. 


Algorithm 1 Pseudo-Optimal Rational Krylov (PORK) 

Input: {E, A, B, C, D), {Sv,R) 

Output: "^2 pseudo-optimal reduced system matrices 
1: V <— AV — EVSy — BR = 0 // Krylov subspace 

2; P~^ = lyap(—S'y, i?^i?) // Low-dimensional Lyapunov eq. 

3; Br = -PrR^ 

4: Ar = Sv+ BrR, Er = /, Cr = CV, Dr = D 


A dual version of the algorithm based on the output Krylov Sylvester equation (4b) 
is available and analogous. The reason why such a reduced model is called H 2 -pseudo- 
optimal is that it is the global optimum, in terms of the "^2 norm, amongst all reduced 
models with the same spectrum. Therefore, the PORK algorithm assigns the eigenvalues 
of the reduced model, resulting from a projection with a Krylov subspace, according to the 
choice of shifts. It is therefore evident that selecting shifts on the right complex half-plane 
yields a stable, pseudo-optimal ROM by construction. However, note that in this setting, 
the choice of appropriate shifts becomes twice as important: not only do they determine 
at which complex frequencies interpolation is achieved, but they also explicitly determine 
the eigenvalues of the reduced model. For this reason, pseudo-optimal reduction has been 
introduced concurrently to an adaptive reduction framework by Panzer et al. in [14], in 
which the choice of shifts for pseudo-optimal reduction is conducted adaptively—through 
a greedy algorithm—within a cumulative reduction framework. A complete and thorough 
exposition of the proposed reduction strategy, known as CUREd SPARK, can be found 
in the recent dissertations of Panzer [11] and Wolf [10] and shall be revised briefly in the 
following. 

Cumulative reduction (CURE) is based on a factorization of the error system that 
allows an iterative reduction of the error term and therefore an adaptive choice of reduced 
order. In fact, if the reduced model, denoted by its transfer function Gr{s), results from 
a projection with an input Krylov subspace V, then the error can be factorized as 


E,A 

B ' 



Br 


■ E,A 

Bi_ 1 


Er, Ar 

Br 

C 

D 

y s 

Cr 

D 

y s 

C 

0 


R 

I 


G(s) Gr(s) Gi(s) GrO 


The first factor G_l(s) represents a high-dimensional system that resembles the original 
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one but is excited only by the input directions not considered during the previous projec¬ 
tion by 11: B± := (/—n)i? = B — EVE~^Bj.. The second factor Gr{s) is low-dimensional 
and resembles the reduced model except for having a unity feedthrough and the output 
matrix corresponding to the matrix of right tangential directions R. 

Using this error factorization, a cumulative reduction procedure is derived which 
reduces, at each iteration k, only the high-dimensional part G±^k of the new error and 
cumulates all ROMs into a total reduced model Gff.: 

G(s) = Gr,l(s) -I- • Grj(s)^ 

Ge,l(s) 

= Gr,l(s) + ^G'r,2(s) + G'_L,2 • Gr,2(s)^ • Gr,l(s) (6) 


— ^r,kM + G_L,fe(s) • G^j,(s). 

A detailed description of the procedure including all expressions, the computation of the 
cumulated system matrices as well as the dual factorization using the output Krylov 
subspace W are given in [11, pp.57-70]. 

Within CURE, there is no restriction of the reduced order at each iteration k. One 
possible choice is to reduce the G_L,fc(s) system at each iteration to a reduced model 
Gr,k{s) of order two. This idea is exploited by the Stability-Preserving Adaptive Rational 
Krylov (SPARK) algorithm, which is aimed at finding a local 772-optimal reduced model 
of order two while reducing the search space to the set of all 'H 2 -pseudo-optimal, stable 
ROMs. Note that this restriction is valid since 'H 2 -pseudo-optimality is a necessary 
condition for 772-optimality, hence the set of all 772-optimal ROMs is entirely contained 
in the set of all 772-pseudo-optimal ROMs. This restriction of the search space has— 
amongst other—the advantage of simplifying the cost function for optimization. In fact, 
if Gr{s) is an 772-pseudo-optimal interpolant of G(s), then the error norm simplifies to 
[10, p.92]: 

||G,(s)||^^ = (G,(s), G,(s))„^ 

= I|G'('S )||«2 - 2 (G(s), Gr(s))^^ -I- ||Gr(s )||^2 (7) 

Since the first term in the resulting expression is constant, the minimization of the error 
norm is equivalent to the maximization of the norm of the reduced model, a cost function 
that is cheap to compute! Restricting the search to 772-pseudo-optimal, stable reduced or¬ 
der models of order two yields simple expressions for gradient and Hessian, parametrized 
by two real numbers. Optimal values can be found by optimization, for example through 
a trust-region descent method [11, pp.75-82]. Cumulation of all 772-optimal reduced 
models of order two yields the final ROM of the desired order. 

2.3. Differential-algebraic equations 

This section revises the main properties of linear, constant coefhcient DAEs and 
introduces the special case of index-1 DAEs in semi-explicit form (SE-DAEs). 
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A DAE as in (1) is primarily characterized by the matrix pair (A,E) defining the 
matrix pencil Pa,eW ■= A + XE. Provided that the pair {A,E) has a finite number 
of generalized eigenvalues, i.e. det Pa,e{X) ^ 0, the pencil is called regular, which is 
often the case if the algebraic constraints are not redundant. In this case the generalized 
eigenvalues can be divided into those that take a finite value and those at infinity. In 
terms of the dynamic system in (1), the former correspond to the dynamic modes whereas 
the latter correspond to the kernel of E, hence the algebraic part. Regularity of the 
pencil also implies that the DAE has a unique solution for a sufficiently smooth input 
and admissible initial conditions, and that there exist two invertible matrices Qi , Qr that 
bring the system to the so called Kronecker- Weierstrafi canonical form [2]: 


Ql E Qr 


Qi A Q. 


0 

N 


i'cyo 



y=[Cf,C^ 

"C^r- 


( 8 ) 


A regular DAE can thereby be decoupled in a dynamical part, represented hy Xf, having 
finite eigenvalues according to the Jordan matrix J, and an algebraic part Xoo with all 
eigenvalues at infinity, described by N, a nilpotent matrix of nilpotency index v. This 
index induces a predominant characterization of DAEs in terms of its “distance” from 
an ODE and is often used to define the type of DAE at hand. Analogously to the 
case of ODEs, asymptotic stability of the DAE can be analyzed by inspecting the finite 
eigenvalues. From the canonical form (8), the transfer function of the DAE can be easily 
computed as [1] 


Gis)=Cf{sInf-J) ^Bf + D + C^ 

Gj{s) >- 



GooO) 


So 


( 9 ) 


It is a characteristic property of DAEs that the transfer function is generally composed of 
a proper part G/(s), resulting from the underlying ODE, and an improper part Goo(s), 
a polynomial in s of degree v — In terms of reduction this has a very important im¬ 
plication: whenever the DAE presents improper behavior, the improper part has to be 
matched exactly to avoid unbounded errors at high frequencies. Generally speaking, the 
correct reduction of a DAE requires the decomposition of the system into dynamic and 
algebraic part, in order to satisfy the algebraic constraints exactly and approximate the 
differential part. However, the computation of the canonical form (8) or, equivalently, 
of projectors onto the respective deflating subspaces, is numerically ill-conditioned and 
not feasible in the large-scale setting. Fortunately, for some DAEs with particular struc¬ 
tures it is possible to identify analytically the relevant subspaces without computing the 
canonical form. 


2.3.1. Semi-explicit, index-1 DAEs 

The special structure of DAE we will consider in this contribution is that of semi¬ 
explicit index-I DAEs (henceforth SE-DAEs), i.e. generalized state-space systems of the 
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form 


Ell 

0 


±1 


All 

Ai2 


Xi 


Bii 

0 

0 


. *2 


A 21 

A 22 


. ^2 _ 

-1- 

B22 


y — [Cii, <^22] 


Xi 

X2 


+ Du 


(10) 


where det^ 22 , detifu 7 ^ 0 , En g and all other matrices partitioned accord¬ 

ingly. Udyn = rankifii is called the dynamic order of the DAE. These systems arise 
frequently in electrical networks and power systems [16-18]. If the semi-explicit form is 
not directly given it can often be achieved by a simple reordering of rows and columns. 
The reason why a DAE as in (10) is called semi-explicit is that the underlying ODE can 
be obtained replacing X 2 through the algebraic equation, which yields: 



Ai Bi 

, -^ ^ 

X\ = All ~ ^ 12 ^ 22^^21 Xi -\- Bii — A 12 A 22 B 22 U 

y = Cii — C22A22 A21 xi + D — C22A22 B22 u. 

■' -V-" '-V-" 

Cl Cl 


( 11 ) 


While the underlying ODE is useful for theoretical considerations, in a large-scale setting 
it is not desirable, if at all feasible, to compute it explicitly. Note that in many applica¬ 
tions, the dimension of the algebraic state-space X 2 is significantly larger than xi. The 
reduction has therefore to be conducted based on the implicit representation ( 10 ). 


2.4- DAE-aware model reduction by Krylov-subspace methods 

In general, the reduction of DAEs cannot be conducted by applying standard methods 
for ODEs. The reason is that while in the ODE case the whole system can be reduced, 
for the DAE it is important to reduce only the dynamic part and keep the algebraic part 
responsible for the improper behavior seen in (9). Therefore, the general approach in 
the reduction of DAEs is to first identify and separate the dynamic from the algebraic 
part and subsequently reduce only the former while preserving the latter. Generally 
speaking, this requires the computation of deflating subspaces and the solution of either 
projected generalized Lyapunov equations (cf. [19]) or the computation of projected 
Krylov subspaces (cf. [20]) depending on the choice of reduction method. A different 
approach that separates algebraic and dynamic part using the concept of tractability 
index is discussed in [ 21 ] and generalized in [ 22 ]. 

This treatise focuses on Krylov-subspace methods due to their generality and compu¬ 
tational advantages. If the DAE presents a special structure, it is generally convenient to 
avoid the explicit computation of the respective subspaces and directly derive algorithms 
that implicitly reduce the ODE. This is done for example by Ahmad and Benner in [23] 
for special index-3 systems and in [24] for special bilinear systems, as well as by Gugercin 
et al. in [20] for special index-I and index-2 systems. Since we will be focusing on SE- 
DAEs, we shall revise in the following the DAE-aware reduction strategy described in 
[20]. This approach is similar to the standard ODE two-sided Rational Krylov and adds 
appropriate shifting terms to ensure interpolation of the underlying ODE and matching 
of the improper part. 
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Theorem 2.1 (Adapted from [20]). Let the full order model he a SE-DAE as in (10). 
Then, the improper part G^o (s) of its transfer function defined in (9) is given by 

Goo{s) = Dimp '■= —C 22 A 22 B22- ( 12 ) 

Eurther, the reduction by tangential interpolation of the underlying ODE can be achieved 
implicitly by computing the reduced model according to 

Er = W^EV, Ar = AV + L Di„,p R 

B, = W^B + LD,mp, Cr = C,V + D,^pR (13) 

Dr — D D^rnp 

where V and W are tangential input and output Krylov subspaces computed with the 
tangential directions in R and L respectively. 

Note that since the DAE considered is of index u = 1, the transfer function cannot be 
improper but can at most contain a constant implicit feedthrough. We shall denote this 
term Dimp to underline that it is the feedthrough implicitly retained in the DAE, as 
opposed to the explicit feedthrough term D. 

We shall revisit this theorem in section 3, where we will give a new proof based 
on a special Sylvester equation. We will also underline the relationship between the 
implicit DAE reduction and the explicit reduction of the underlying ODE, as this is of 
central importance for the fidelity of the ROM. In this context, it will be shown that the 
reduction proposed in theorem 2.1 cannot be applied for orthogonal projections V=W 
of arbitrary SE-DAEs. 

2.5. Problem statement 

The goal of this contribution is to introduce techniques that can be used to ensure 
the stability of the ROM obtained by Krylov-based reduction of SE-DAEs. This is done 
by extending the concepts and methods introduced in this section to this class of DAEs. 
Even though the structure of SE-DAEs suggests how to compute the underlying ODE 
explicitly, the computations required for the DAE-aware, stability preserving reduction 
will use the original DAE matrices only. It will be shown how the proposed procedures 
are equivalent to a the corresponding reduction performed directly on the underlying 
ODE. Numerical investigations on different benchmark systems will be used to assess 
the effectiveness of the proposed methods. 

3. Stability-preserving reduction for strictly dissipative DAEs 

The first case of stability preserving reduction is suitable for SE-DAEs that are given 
in a strictly dissipative realization. Similarly to the ODE case, orthogonal projection will 
be used to preserve dissipativity. Nevertheless, note that orthogonal projection cannot 
be conducted correctly on any SE-DAE, even by applying the DAE-aware MOR scheme 
of theorem 2.1. In general, the resulting ROM will fail to approximate only the dynamic 
part and lose dissipativity. The conditions the DAE has to satisfy in order to obtain a 
correct reduction in this sense will be derived in the following. 
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3.1. Strictly dissipative SE-DAEs 

The concept of a strictly dissipative formulation introduced in section 2.2.1 for ODEs 
can be naturally extended to SE-DAEs by inspection of the underlying ODE, as stated 
in the following definition 

Definition 1. A SE-DAE System as in (10) with underlying ODE (11) is said to be in 
strictly dissipative form if and only if 

El = E^ 0 A Ai A^ A 0. 

Clearly this definition respects the same properties as the one in section 2.2.1 and, in 
particular, it implies stability of the SE-DAE. Note that in some cases strict dissipativity 
of the SE-DAE can be assessed without explicitly computing the underlying ODE, as it 
is explained in the next proposition. 

Proposition 1. Let A and Ai be defined as in (10) and (11) respectively. Assume 
A + A^ -< 0. Then Ai + Aj -< 0. 

Proof. The proof is given in Appendix A. □ 


Obviously, not all SE-DAEs are given in strictly dissipative form. Nonetheless, if the 
SE-DAE is asymptotically stable, then analogously to the ODE case, a transformation 
to strictly dissipative form can be found: 

Lemma 1. Assume the SE-DAE in (10) is asymptotically stable. A strictly dissipative 
formulation is obtained by multiplying the DAE from the left with the matrix 


EjiP -EjiPAi^Af^ 
0 I 


where P= aO solves the Lyapunov inequality 

EjiPAiP AjPEii A 0 

and Ai is the Schur complement defined in (11). 

Proof. Note that the transformation T can be factorized as 


EjiP 

1 

o 


I —^12^22^ 

1 - 

o 

1 


0 / 


Applying this transformation to the SE-DAE yields 


EjiPEii 

0 

X = 

Ej^PAi 0 

X + 

EjiPBi 

0 

0 


A 21 A 22 


P 22 


which is in strictly dissipative form by construction. □ 


Therefore, finding a strictly dissipative transformation reduces to solving a convex fea¬ 
sibility problem, which might still be difficult in the large-scale case. Alternatively, a 
Lyapunov equation with appropriate right-hand side might be solved. Note at this point 
that if the FOM is not given in strictly dissipative form, stability preserving reduction 
can still be achieved by using the second method proposed in section 4. 
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3.2. Orthogonal projection of SE-DAEs 

Combining what has been discussed so far, we are now ready to establish the model 
reduction strategy that acts only on the sparse matrices of the SE-DAE and ensures 
stability of the ROM. Similarly to the ODE case, we will approach the problem by 
applying orthogonal projection to the DAE. First, we start by demonstrating how the 
reduction framework as in theorem 2.1 generally fails to comprise the special case of 
orthogonal projection. This is counterintuitive, since for ODEs orthogonal projection is 
a special case of skew projection (V W). To do so, we compare from a theoretical 
standpoint the reduction of the underlying ODE to the implicit reduction of the DAE. 

Recalling the SE-DAE (10) and the corresponding explicit representation of the un¬ 
derlying ODE (11) we start by showing the equivalence of the Krylov subspaces computed 
with either representation. 

Lemma 2. Given the SE-DAE (10) and the underlying ODE (11), the following equiv¬ 
alence between Sylvester equations holds 


A 


Vi 

V2 


- E 


Vi 

V2 


Sv-BR = 0 


AiVi - EiViSv - BiR = 0 


(14a) 

(14b) 


where V = [V^, ^ 


is partitioned according to Udyn- 


Proof. 


A 


^2 



Sv-BR = 0 


f AnVi-h A12V2 - EnViSv - BnR = 0 
( V2 = A22 (—A21V1 -|- B22R) 


AiVi 


EiViSv - BiR = 0 


□ 


Naturally, the dual result in terms of the output Sylvester equation (4b) holds as well. 
This result implies that the projection matrices Vi and Wi needed to approximate the 
underlying ODE can be computed by solving large sparse system of equations in terms of 
the original DAE matrices without explicitly computing Schur complements. Even more 
importantly, this result implies, through the duality between Krylov and Sylvester, that 
both the reduced model obtained by reducing the underlying ODE using Vi, Wi and the 
one obtained through the reduction of the DAE using V, W share the same interpolation 
data {Sy, R) and hence achieve tangential matching of the same moments of the original 
model. However, this result does not guarantee that these reduced models will actually 
be the same. In order to get this result, we have to ensure that the second procedure 
implicitly reduces the underlying ODE by operating on the DAE. We can analyze this by 
projecting the ODE and SE-DAE with their respective projection matrices and compare 
the expressions. The results are summarized in table 1, where the last column includes 
the correction terms of theorem 2.1. 

The results in the table show that the two reduced models differ exactly by the 
terms that are compensated in the SE-DAE aware reduction algorithm of theorem 2.1, 
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ODE {Vi,Wi) 

SE-DAE (E, W) 

correction [20] 

Er 

WffEiVi 

WffAiVi 

WffEiVi 

LDimpR 

Br 

W^Bi 

Bi —LDimp 

LDimp 

Cr 

CiVi 

Cl El —DimpR 

Dimp LI 

Dr 

B 4“ E)irnp 

D 

Dimp 


Table 1: Comparison of a skew projection for ODE and SE-DAE 


confirming that this procedure effectively yields moment matching and implicit reduction 
of the underlying ODE for a skew projection. 

Will this result hold true also for one-sided reduction? The theory on ODE reduction 
suggests that this should be the case. However, it turns out that this is not always true 
when operating with SE-DAEs, as we will state in the following theorem. 

Theorem 3.1. Consider a SE-DAE as in (10) and its underlying ODE defined in (11). 
The ROM obtained by applying orthogonal projection (W = V) of the SE-DAE, being 
V a basis of the Krylov subspace as in (14a), is equivalent to the one resulting from the 
respective orthogonal projection applied to the underlying ODE with Vi satisfying (14b), 
provided that B 22 = 0. 

Proof. The proof is obtained by straightforward computation of the orthogonal pro¬ 
jection of the SE-DAE with V = s-nd by using the relationship V 2 = 

Af 2 (— A 21 V 1 -I- B 22 R) obtained in lemma 2. The correction terms of table 1 will be 
applied. Since we only compute V as an input Krylov subspace, the matrix of left tan¬ 
gential directions L is set to 0. Subsequently, the ROM obtained is compared to the one 
that would result from a direct orthogonal projection of the underlying ODE. For brevity 
we will omit the computations and directly compare the reduced matrices in table 2. 



ODE (lEi=Ei) 

SE-DAE (1E=E), corrected [20] 

Er 

V^EnVi 

C/CiiEi 

Ar 

V^AiVi 

C/ Ai El 

Br 

V^Bi 

ky (^11 — ^ 21 ^ 22 ^ B 22 ) + R^ BJ 2 A 22 B 22 

Cr 

Cl El 

Cl El 

Dr 

B 4“ Birap 

D DfjYip 

Table 2: Comparison of an orthogonal projection for ODE and SE-DAE 

The results show how this procedure generally fails to correctly reduce the underlying 


ODE. Nevertheless, if the original model satisfies B 22 = 0, then the reduction is correct 
after all. □ 


Accordingly, care needs to be taken when applying orthogonal projection to reduce 
DAEs. Even though moment matching is still guaranteed by the computation of the 
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Krylov subspace, the reduction fails to project the underlying ODE unless the condition 
B 22 = 0 is satisfied. Even more importantly: the additional terms in table 2 might cause 
loss of dissipativity of the reduced model! The relevance of this will be shown through 
a numerical example in section 5. Further note that the assumption B 22 = 0 is relevant 
from a practical standpoint but restrictive from a theoretical one, since the results in 
table 2 infer that a system satisfying A 22 = AJ 2 , ^12 = Aj^, C 22 = B 22 could be correctly 
reduced through orthogonal projection by choosing L = BA in the correction terms. In 
fact, in this case BJ 2 A 22 B 22 = —LDimp and Bn — AI 1 A 22 B 22 = Bi. 

Finally, note that this result extends naturally to the dual case of output-based or¬ 
thogonal projection. 

Corollary 1. Consider an SE-DAE as in (10) and its underlying ODE defined in (11). 
The ROM obtained by applying orthogonal projeetion (V = W) of the SE-DAE, being W 
a basis for the Krylov subspace as in (4b), is equivalent to the one resulting from the 
respective orthogonal projection applied to the underlying ODE, provided that C 22 =0 or 
A 22 = AJ 2 ,Ai 2 = AJ^ , C 22 = BJ 2 ■ 

Proof. The proof is analogous to the dual case. Note that also in this case, the condition 
C 22 = 0 can be dropped if the system satisfies A 22 = AJ 2 , A 12 = Af ^, C 22 = BJ 2 and the 
right tangential directions for the correction terms are chosen such that R = LA . □ 

To sum up, provided that the SE-DAE system considered is not excited through 
the algebraic equations or the output influenced by algebraic variables, then DAE-aware 
reduction by orthogonal projection as proposed in theorem 2.1 can effectively reduce 
the underlying ODE implicitly, by choosing the right type of Krylov subspace (input 
or output). If the system has some symmetries with respect to the algebraic variables, 
namely A 22 = AJ 2 , A 12 = Aj^, C 22 = BJ 2 , then both input and output Krylov subspaces 
can be used, provided the tangential directions are chosen such that R = LA . 

3.3. Stability-preserving reduction by orthogonal projection 

Due to the previous results, we understand that we can implicitly apply orthogonal 
reduction to the underlying ODE choosing either input or output Krylov subspaces de¬ 
pending on the structure of the system. As we know from ODE theory, this process 
preserves stability in case that the underlying ODE—and hence by definition 1 also the 
SE-DAE—is given in strictly dissipative form. 

Theorem 3.2. Assume the SE-DAE as in (10) is given in strictly dissipative form. 
Further assume that either B 22 = 0 or A 22 = AJ 2 , A 12 = Aji, C 22 = Then reduction 

of the SE-DAE through orthogonal projection using the input Krylov subspace defined in 
(14a) and the DAE-aware procedure of theorem 2.1 yields an asymptotically stable ROM 
that implicitly reduces the underlying ODE. 

Proof. The proof is obtained by observing that due to lemma 2, the proposed reduction 
of the SE-DAE is equivalent to the direct reduction of the underlying ODE. Applying 
orthogonal projection preserves strict dissipativity and hereby stability. □ 

Naturally the dual result holds as well. 
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Corollary 2. Assume the SE-DAE as in (10) is given in strictly dissipative form. Eur- 
ther assume that either C 22 = 0 or A 22 = ^^ 21^12 = ^ 21^^22 = - 822 - Then reduction 
of the SE-DAE through orthogonal projection using an output Krylov subspace and the 
DAE-aware procedure of theorem 2.1 yields an asymptotically stable ROM that implicitly 
reduces the underlying ODE. 

The effectiveness of this method will be assessed through numerical examples in sec¬ 
tion 5. In addition, it will be shown that if the correct choice of Krylov subspace is 
ignored, the ROM might indeed become unstable. 

4. Stability preserving, adaptive reduction by 'H 2 -pseudo-optimality 

The preservation of stability in the case of ODEs can be achieved by construction 
using 'H 2 -pseudo-optimal reduction as presented in section 2.2.2. In this procedure, if the 
shifts of the Krylov subspaces are chosen on the right complex half-plane, then the ROM 
will have all eigenvalues in the left complex half-plane and hence be asymptotically stable. 
In the following, we will extend this result to the case of SE-DAEs. Since the choice of 
shifts becomes twice as important in the pseudo-optimal setting, we will subsequently 
address the question of how to appropriately select the shifts. 

4-.1. 'H 2 -pseudo-optimal reduction of SE-DAEs 

Recall the construction of pseudo-optimal ROMs following the PORK algorithm 
(cf. algorithm 1). One peculiar characteristic of PORK is that the reduced matrices 
{Er, Ar, Br) —Or {Er, Ar,Cr) in the dual version—are independent of the original model 
and depend only on the interpolation data, i.e. merely on the reduced eigenvalues and 
tangential directions encoded in the pair {Sv,R) or {Sw,L) respectively. This char¬ 
acteristic makes it particularly suitable for the extension to SE-DAEs, at least in case 
the reduced model should be an ODE and not a DAE. Note however, that since the 
DAEs considered are of index v = 1, their transfer behavior can be modeled entirely by 
an equivalent ODE, as it can be seen from the transfer function in (9). Therefore, the 
restriction to reduced ODE models can be conducted without loss of generality. 

It turns out that the PORK algorithm can be adapted to SE-DAEs by combining 
it with the appropriate correction terms proposed in theorem 2.1. Since the proof is 
straightforward following what has been said so far, we will limit the exposition to the 
adapted algorithm. 


Algorithm 3 SE-DAE PORK 
Input: {E, A, B, C, D), {Sv,R) 

Output: 'H 2 -pseudo-optimal reduced system matrices 
1: V ^ AV — EVSv — BR = 0 // Krylov subspace 

2 ; P~^ = lyap(—// As in the ODE case 
3: Br = -PrR^ 

4: Ar = Sv + B^R, Ey = I 

5 : Cr = CV + DimpR, Dj. = D -\- Dimp j j cf. theorem 2.1 and (12) 


Naturally this result extends to the dual case. This stability preserving reduction 
procedure clearly is far less restrictive than the one proposed in section 3. On the one 
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hand, it does not require the original model to be given in strictly dissipative form. On 
the other hand, it poses no restrictions on the structure of the original model and can 
therefore be applied to any SE-DAEs. 

What still remains unclear at this point is how to choose appropriate shifts. Recall 
that in standard reduction based on Krylov subspaces the choice of shifts only affects 
the matching frequencies, while the eigenvalues implicitly result from the reduction. 
By contrast, in the pseudo-optimal setting the reduced eigenvalues are directly chosen 
through the choice of shifts. While it is known that pseudo-optimality is a necessary 
condition for H 2 -optimality (cf. Maier-Luenberger conditions [8]), this does not imply 
that pseudo-optimal ROMs are better than the equivalent ROM achieved, for instance, 
through two-sided reduction. The choice of shifts becomes even more crucial and the 
next section will show how these can be determined adaptively even for SE-DAEs within 
the reduction framework CUREd SPARK. 

J^.2. Judicious choice of shifts: CUREd SPARK for SE-DAEs 

Analogously to the ODE case, the cumulative reduction framework with adaptive 
choice of shifts (CUREd SPARK) can be extended to the case of SE-DAEs to complement 
the pseudo-optimal reduction with a judicious choice of shifts. With this respect, it is 
worth separating the discussion for the two complementary classes of SE-DAEs: the ones 
presenting and the ones omitting the implicit feedthrough term D^^p. 

Whenever the SE-DAE does not have an implicit feedthrough term Dimp as defined 
in (12), then the original CUREd SPARK introduced in section 2.2.2 can be used without 
restrictions. In fact, the algorithms do not require the matrix E to be regular. Numerical 
examples in section 5 will show the effectiveness of this procedure. It is worth noting at 
this point that the vast majority of benchmark SE-DAEs models commonly available fall 
into this category. This is motivated by physical intuition, since real technical systems 
mostly have some delay between excitation and response and are hence feedthrough-free. 

On the other hand, if an implicit feedthrough term Dimp is present, then an adapta¬ 
tion of CUREd SPARK is required. It can be shown that the iteration in CURE (cf. (6)) 
can be adapted to take care of the implicit feedthrough term Dimp resulting at each it¬ 
eration k. However, the greater issue arises when applying SPARK to SE-DAEs and 
requires a modification of the DAE, as we shall discuss in the following. 

As seen in section 2.2.2, the goal of SPARK is to find a reduced order model Gr(s) 
of order n=2 that is locally a 'H 2 -optimal approximation of the original model G(s). 
Exploiting pseudo-optimality, the error norm is given by (cf. (7)) 

||Ge(s)||^/=="°||G(s)||^^-||G.(s)||^^ 

This was true at least in the ODE case, where it is common practice to leave feedthrough 
terms D aside from the reduction and integrating them subsequently into the reduced 
model. However, in the case of DAEs retaining an implicit feedthrough term Dimp, this 
term is hidden inside the matrices A, B, G and cannot be directly removed prior to the 
reduction. In this case, obviously the above error expression becomes meaningless since 
the 'H 2 -norm of a non strictly proper system does not exist. To better understand this 
fact, let us separate the strictly proper part G®p(s) of the transfer function G(s) from 
the implicit feedthrough Dimp and assume that the reduced model Gr{s) is also given 
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as the sum of a strictly proper part G^^(s), which is a pseudo-optimal approximant of 
G^P{s), and the same feedthrough term Di^np- The error expression then becomes 

\\Ge{s)\\l^ = {G{s)-Gr{s),G{s)-Gr{s))^^ 

= {G^^{s) + Amp - GlP{s) - Amp, G(s) - G,(s))„^ 

\\G^^{s)\\l,^-\\G7’{s)\\l^ (15) 

In this case, the feedthrough term cancels out and by pseudo-optimality we obtain the 
difference of two well defined norms. Therefore, in order for (15) to hold true, we need 
to: 


a) find an expression for G^^(s), the strictly proper part of G(s), 

b) find a pseudo-optimal approximation G®^(s) for G'*p(s), 

c) make sure the reduced model Gr{s) retains the feedthrough term Dimp, i.e. Gr{s) = 
GfF{s) + Dimp- 


Whenever the SE-DAE has an implicit feedthrough term Dimp, the model reduc¬ 
tion procedure has to compute it and take it into consideration during reduction (cf. 
theorem 2.1 or [20]). Therefore, its inclusion in the reduced model is easy to implement. 
Further, provided a strictly proper model G®p(s) is given, all the algorithms discussed so 
far can be implemented to reduce it. It follows that the only open question at this point 
remains how to extract a strictly proper representation G^^(s) out of an SE-DAE with 
implicit feedthrough. One possible solution is presented in the following: 

Proposition 2. Consider a SE-DAE as in (10) with an implicit feedthrough term Dimpf^ 
0, whose transfer function can he written as the G{s)=G'^P{s)-\-Dimp- Then a state-space 
realization of the strictly proper part G^p{s) is given by 


Ai 

0 ■ 


±1 


■ An 

Ai2 


Xi 

_1_ 

Bll — ^12^22^ B22 

0 

0 


X 2 


A 21 

A 22 


. ^2 


0 


y — [Gii, G 22 ] 


(16) 


or alternatively 


Ell 

0 ■ 


±1 


■ An 

Ai2 


Xi 

0 

0 


. *2 


A21 

A22 


X2 


y — [Gii — G22A22*‘A21,0] 


Ai 

A 2 


u 


Xi 

X2 


(17) 


Proof. The proof is straightforward and follows by computing the underlying ODE for 
both the system in (16) and (17) as it was done in (11). □ 

Note that, if an implicit feedthrough term Amp = —G 22 Af 2 B 22 is present in the 
system, then the reduced model must pertain it in order to have a bounded approximation 
error. This implies the fact that the term Dimp must be computed anyways. Therefore, 
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the vectors ^^ 2^,622 or C 22 A 22 required for computing the strictly proper realization 
G^P{s) are already available without further computations. 

Exploiting this result, it is possible to choose reduced order and interpolation points 
adaptively by means of CUREd SPARK even for SE-DAEs with implicit feedthrough 
term Dimp- We will illustrate the effectiveness of this approach through numerical ex¬ 
amples in the following section. 

5. Validation through numerical examples 

The theoretical results presented in the previous sections are validated through some 
numerical examples. First, we introduce the models used for the computations, then we 
discuss the results. 

5.1. Transmission line model 

A transmission line, i.e. the transmission of electrical signals through a one-dimensional 
conductor, is described by a partial differential equation (PDE) in both time and space. 
A common way of approximating this PDE is by using basic elements of an electric net¬ 
work distributed along the line, discretizing the PDE in space [25]. A typical, simplified^ 
representation of such a discretized model is given in hgure 1 . 



Figure 1: Electrical circuit of a transmission line approximation 

The input of the system is represented by a voltage source C/q. represent 

the distributed resistance, inductance and capacitance of the line respectively. The typ¬ 
ical output of the system is the voltage over the capacitor at the end of the line. 
The number of loops used to discretize the transmission line shall be denoted by q. In 
general, the higher the number of loops, the better the approximation. The values for 
the distributed parameters are taken from [25, p.588]. In general, the parameters vary 
largely depending on the geometry and materials of the problems. For our simulations, 
the special case of a telephone line and a transmission frequency of 1 Hz was taken, for 
which the parameters are given in table 3. 

The modeling is conducted using hrst principles, such as Kirchhoff’s laws and consti¬ 
tutive equations for the elements (cf. e.g. [26]). Each loop is modeled using the variables 
Xi = [iRi, Uci,UR., Ici, ULi]^, which can be stacked up to yield the Sq-dimensional state 


^The simplification results form the omission of the conductance that is usually shunt in parallel 
to the conductances in such a network and is justified by the selection of parameters (cf. [25, p.588]) 
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Parameter 

Value 

R'r 

172.24 • 10-3 n/m 

L) 

0.61 • 10"® 

ci 

51.57- 10-32 F/^ 


Table 3: Transmission line parameters used for simulations 


vector that characterizes the DAE. The general structure of the SE-DAE matrices re¬ 
sulting from this modelling can be found in Appendix B. Note that for the standard 
choice of input and output as in figure 1 the system matrices satisfy B 22 7 ^ 0 (the input 
Uo enters the equations through Kirchhoff’s law, an algebraic constraint), (722 = 0 (the 
output is a dynamic variable) and A 22 7 ^ ^^ 2 - Also note that even though the modeling 
does not yield directly a strictly dissipative representation, the procedure described in 
theorem 1 was used to find a suitable representation up until an order of A^=700, which 
corresponds to <7=140 loops. Einally, note that by appropriately selecting the output of 
the system, e.g. by choosing the voltage over the first inductance L'^, it is also possible 
to model a system with implicit feedthrough Dimp- 

5.2. Power system benchmark models 

Another set of benchmark systems chosen to test the validity of the proposed algo¬ 
rithms is given by the power system examples created at the Brazilian Electrical Energy 
Research Center (CEPEL) and available online in the MOR Wiki [18]. These systems 
represent large power systems (including lines, buses, power plants etc.) linearized about 
an operating point and are used to simulate and study the oscillations of complex power 
systems. The reduction of such systems is relevant for the numerical simulations required, 
for instance, for small-signal stability analysis, controller design, and real-time investi¬ 
gations of the transient behavior. Eor more informations on the origin of the systems, 
please refer to [17, 18]. 

Most of these systems are index 1 DAEs that are either already in semi-explicit form 
or can be transformed to it by simple reordering of rows and columns. The models 
are mainly strictly proper, meaning that they have neither explicit feedthrough D nor 
implicit feedthrough Dimp- However, the so called “MIM046” system has SISO transfer 
functions on the diagonal of the transfer function matrix that present such an implicit 
feedthrough 

For the numerical investigations of this treatise, two different models of the Brazil¬ 
ian Interconnected Power System of 1997 (BIPS/97) were used. Both of them have 
a state-space dimension of roughly N=13250 and a dynamic order, that is the order 
of the underlying ODE, of ndy„=1664. While the first one is SISO and strictly proper 
{wwjuref JSA05.mat), the second one is MIMO with implicit feedthrough on the diagonal 
channels {mimoAQxA(S_system.mat). 

5.3. Results 

5.3.1. The importance of SE-DAE-awareness in orthogonal reduction 

The first example is aimed at showing what can happen if the results of section 3.2 
are neglected. Figure 2 shows the results obtained on a transmission line model with 
<7 = 10 loops. 
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Figure 2: Orthogonal projection using input and output Krylov subspaces for the transmission line 
example {q = 10) 

The bode plot of the FOM with N=50 is compared to two ROMs of order n=20, one 
obtained through V-based and the other through W-based orthogonal projection of the 
SE-DAE. Recall that since C 22 =0, B 22 y^O and A 22 ^ AJ 2 , only W-based orthogonal 
projection is guaranteed to effectively reduce the underlying ODE. The shifts for the 
Krylov subspaces were chosen along the imaginary axis with frequencies corresponding 
to the peaks of the FOM. As it can be seen, the W-based ROM matches the FOM 
perfectly. In fact, 20 corresponds to the dynamic order of the system, i.e. the order of 
the underlying ODE. Therefore, the W-based ROM is merely a transformation of the 
underlying ODE and hence shares the same transfer function. On the other hand, the 
figure shows how this is not true for the V-based ROM. Moment matching at the shifts 
still holds, however, the reduction results are clearly worse. In particular, this is not a 
transformation of the underlying ODE anymore. 

5.3.2. Stability preserving reduction of strictly dissipative SE-DAEs 

The effectiveness in preserving stability while interpolating the underlying ODE of 
the reduction method proposed in theorem 3.1 is shown on a transmission line model 
with <7=140 loops and state-space dimension V=700. Before reduction, the SE-DAE has 
been transformed to a strictly dissipative formulation as explained in theorem 1. The 
results are given in figure 3. The FOM is compared with two ROMs of order n=100, one 
obtained with V-based, the other with W-based orthogonal projection about the origin. 
Recall that this system satishes C 22 = 0, hence only IV-based orthogonal projection is 
expected to yield acceptable results. 

The plot suggests that the W-based ROM yields a slightly better approximation of the 
FOM at higher frequency range. The plot does not show that while the W-based ROM 
coincides with an equivalent ROM obtained reducing directly the underlying ODE, the V- 
based one does not, emphasizing again the importance of SE-DAE-aware reduction. Even 
more importantly: while the W-based ROM preserves strict dissipativity and stability, 
the V-based ROM is unstable! These results are summarized in table 4. 
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Figure 3: Orthogonal reduction using input and output Krylov subspaces for the transmission line 
example {q = 140) 


dissipative 

stable 


FOM ROM (W) ROM (V) 

“7 7 X 

/ / X 


Table 4: Comparison of different orthogonal projections of the transmission line 


5.3.3. Stability preserving, l-L 2 -pseudo-optimal reduction of general SE-DAE 

The stability preserving, adaptive reduction of the SE-DAE with CUREd SPARK 
is possible even for systems that are not in strictly dissipative form. As discussed in 
section 4, whenever the system has no implicit feedthrough Dimp, CUREd SPARK can 
be applied by simply replacing the PORK algorithm for 'H 2 -pseudo-optimal reduction 
by the SE-DAE PORK (algorithm 3). The reduction results using this strategy for the 
strictly proper power system described in 5.2 are shown in figure 4, where the magnitude 
plot of the FOM (in dB) is compared to the magnitudes of the error systems resulting 
from a. W OT U-based reduction through CUREd SPARK respectively. 



Figure 4: CUREd SPARK reduction of the BIPS/97 system (Dimp = 0, n = 50) 
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The original order of the DAE is A^=13251, the dynamic order is ridyn = 1664. The 
reduction was conducted using CUREd SPARK to a reduced order of n=50, both with 
V and IE-based factorization. The selection of the shifts was conducted automatically 
in each iteration of CURE by the greedy, 'H 2 -optimal algorithm in SPARK. As it can 
be seen from the magnitude of the frequency response, both reduced models achieve a 
satisfactory approximation of the POM while preserving stability. 

The results of the reduction of SE-DAEs presenting an implicit feedthrough are shown 
in figure 5 and figure 6. The FOM represents one of the diagonal input-output channels of 
the MIMO power system described in section 5.2. Like the previous system, the original 
order is iV=13250 and the dynamic order is ndy„=1664. The reduced order was chosen 
to n=50 for all channels considered. Also in this case, the choice of shifts was conducted 
automatically at each iteration of CURE by the greedy algorithm in SPARK. 
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Figure 5: CUREd SPARK reduction of the MIMO BIPS/97, channel 25 (Dimp 0,n = 50) 
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Figure 6: CUREd SPARK reduction of the MIMO BIPS/97, channel 42 (Dimp ^ 0, n = 50) 


The figure shows how the proposed reduction algorithm manages to detect the relevant 
dynamics of the systems and preserve both stability and the implicit feedthrough. 
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6. Conclusions 


In this contribution the stability preserving reduction of differential-algebraic equa¬ 
tions (DAE) of index 1 in semi-explicit form has been discussed. To the authors’ knowl¬ 
edge, this is the first work to address this issue. To this purpose, two different methods 
have been presented. 

The first is aimed at preserving the strictly dissipative form of the underlying or¬ 
dinary differential equation (ODE) by means of orthogonal projection of the DAE. A 
transformation of the DAE to strictly dissipative form through the solution of a convex 
feasibility problem has been introduced. It has been shown in theory and through nu¬ 
merical examples that the correct orthogonal projection can be achieved only if the DAE 
fulfills certain conditions. In this case, the appropriate choice of input or output Krylov 
subspace for the orthogonal projection yields a strictly dissipative, hence stable reduced 
order model. 

The second method uses 'H 2 -pseudo-optimal reduction and is more general in that 
it does not impose any conditions on the DAE. For this purpose, the Pseudo-Optimal 
Rational Krylov (PORK) algorithm was adapted to include this class of DAEs. By 
construction, H 2 -pseudo-optimal reduction generates asymptotically stable reduced order 
models if the shifts of the Krylov subspaces are chosen on the right complex half-plane. 
However, in this setting the choice of appropriate shifts is even more important for a 
good approximation quality. If the DAE has no implicit feedthrough term, which is 
often the case for technical systems, then this adapted PORK algorithm can be used 
within CUREd SPARK to adaptively select reduced order and shifts. If the system at 
hand does possess an implicit feedthrough term, then a realization of the strictly proper 
part was presented. This strictly proper part can be reduced adaptively and stability 
preserving with CUREd SPARK and the complete feedthrough term is added at the end. 

The proposed methods have been tested numerically showing their validity on differ¬ 
ent academic and real-life systems. 


Appendix A. Proof of proposition 1 


Proof. The proof of proposition 1 amounts to showing that the Schur complement of a 
strictly dissipative matrix, i.e. satisfying A+A^ -< 0, is strictly dissipative itself. For this 
purpose, we shall partition the matrix A according to (10) 


A = 


All Ai2 
A21 A22 


and define the Schur complement Ai := An — Ai 2 A^ 2 ^A 2 i. It then follows 

A -I- A^ ^ 0 A~^ (A -I- A^) A“^ = A“^ -|- A~^ -< 0 

where the invertibility of A follows from its strict dissipativity. The inverse of A can be 
build block-wise as follows (* denotes entries that are not required for our reasoning) 


A-i = 


AI 

* 


-1 


(A.l) 
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Appendix B. Transmission line SE-DAE 

In the following we introduce the model of the transmission line presented in figure 
1. Since the aim is to analyze the characteristics of the DAE-aware procedures presented 
in this work, the model has been kept as simple as possible. The conductance Gi that is 
typically in parallel to the capacitance Ci has been neglected. This case corresponds e.g. 
to the transmission line model of a 24 gauge telephone PIC cable at 21°C' and a 1 Hz 
frequency [25]. The equations have been derived using first principles of physics due to 
their simplicity and direct physical interpretation [26]. Each loop i out of q has a set of 
5 states, namely Xi = [lRi,Uci,Uii-, Ici,ULi]^, where we have already made use of the 
fact that the current flowing through the resistances and inductances in each loop must 
be identical. The five equations that model each loop and define the relations amongst 
the states are Kirchhoff’s laws and the constitutive equations of each element (resistor, 
inductor, capacitor). For each loop i = 2,... ,q it therefore holds 


+ Uc, - Uc^_, = 0 

(B.la) 

O 

II 

1 

+ 

1 

(B.lb) 

C/fl. - ITJr, = 0 

(B.lc) 


(B.ld) 


(B.le) 
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For z = 1, equation (B.la) is changed replacing C/ci_i by the input C/q- This formalism 
leads to the SE-DAE defined block-wise as 
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(B.2) 


where Iq is the identity matrix of size q, Om,n is a zero matrix of size mxn, Iq and /“ 
are squared matrices having ones on the first super- and sub-diagonal respectively and 
the other entries are defined as follows 

Lq := diag(Li,...,Lg) 

Ir '■= [Iri^ ■ ■ ■ 

:= [1,0,... 0]^ 

(5c := [0,...,0,1]. 


Note that other techniques like Modified Nodal Analysis (MNA) are widely used to 
model such networks [2]. However, they do not directly yield a SE-DAE as in (10) and 
were therefore deemed as less suitable for the purpose of this treatise. It is however 
possible to transform the resulting equations from said modeling formalisms into a SE- 
DAE, as the underlying system stays the same. 
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